% Table 3 in "Commodity Prices, Convenience Yields, and Inflation,"
% March 2011

clear all;

load cyj.dat;                   % convenience yields on 23 commodities
load us_intrate_monthly.dat;    % US 3-month T-bill rate
load us_cpi_monthly.dat;        % US CPI indeces (all items, lfe etc.)
load imf_index_monthly.dat;     % IMF aggregate commodity price index

group1=4; group2=10; group3=12; group4=16; group5=21; group6=23;
i=us_intrate_monthly;
P=us_cpi_monthly;
SA=imf_index_monthly;
nn=rows(P);

inflmeas=1;     % inflation measure (all items, lfe etc.)
rhat=2;         % number of principal components
nwl=-1;         % number of lags in Newey-West HAC estimator
bsize=4;        % bootstrap block size
B=4999;         % number of bootstrap replications

[ehat,Fhat,lamhat,ve]=pc(standard(cyj),rhat);  % PC estimation

hh=[1 3 6 12]';  % forecast horizon

for hi=1:rows(hh);
h=hh(hi);

infl_h=100*(log(P(3+h:nn,inflmeas))-log(P(3:nn-h,inflmeas)));
infl_1=100*(log(P(3:nn-h,inflmeas))-log(P(2:nn-h-1,inflmeas)));
infl_2=100*(log(P(2:nn-h-1,inflmeas))-log(P(1:nn-h-2,inflmeas)));
p1_cy=Fhat(3:nn-h,1);
p2_cy=Fhat(3:nn-h,2);
cy_ave1=mean(cyj(3:nn-h,1:group1)')';
cy_ave2=mean(cyj(3:nn-h,group1+1:group2)')';
cy_ave3=mean(cyj(3:nn-h,group2+1:group3)')';
cy_ave4=mean(cyj(3:nn-h,group3+1:group4)')';
cy_ave5=mean(cyj(3:nn-h,group4+1:group5)')';
cy_ave6=mean(cyj(3:nn-h,group5+1:group6)')';
ds_1=100*(log(SA(3:nn-h,1))-log(SA(2:nn-h-1,1)));

% Baseline inflation equation
x=[ones(rows(p1_cy),1) p1_cy p2_cy cy_ave1 cy_ave2 cy_ave3 cy_ave4 cy_ave5 cy_ave6 ds_1 infl_1 infl_2];
res=nwest(infl_h,x,0);
b=res.beta;
r=res.resid;
R2=res.rbar;
G=nw((r*ones(1,cols(x))).*x,nwl);
V=inv(x'*x)*G*inv(x'*x);
se = sqrt(diag(V));
ts = b./se;

x2=[ones(rows(p1_cy),1) ds_1 infl_1 infl_2];
res=nwest(infl_h,x2,0);
b2=res.beta;
R2_2=res.rbar;

% Bootstrap procedure
bmatrix=[cyj(3:nn-h,:) infl_h infl_1 infl_2 p1_cy p2_cy ds_1];
gg=floor(rows(bmatrix)./bsize);
bmatrix=bmatrix(1:gg*bsize,:);
betab=[]; betab1=[]; betab2=[];
tstatb=[]; tstatb1=[]; tstatb2=[];
for j=1:B;
    indexx=floor(1+gg*unifrnd(0,1,1,gg))';
    ind=blockind(indexx,bsize);
    ind=ind(1:rows(bmatrix));
    matrixb=bmatrix(ind,:);
    cyjb=matrixb(:,1:23);
    inflb_h=matrixb(:,24);
    inflb_1=matrixb(:,25);
    inflb_2=matrixb(:,26);
    p1_cyt=matrixb(:,27);
    p2_cyt=matrixb(:,28);
    dsb_1=matrixb(:,29);
    cyb_ave1=mean(cyjb(:,1:group1)')';
    cyb_ave2=mean(cyjb(:,group1+1:group2)')';
    cyb_ave3=mean(cyjb(:,group2+1:group3)')';
    cyb_ave4=mean(cyjb(:,group3+1:group4)')';
    cyb_ave5=mean(cyjb(:,group4+1:group5)')';
    cyb_ave6=mean(cyjb(:,group5+1:group6)')';
    
    [ehatb,Fhatb,lamhatb,veb]=pc(standard(cyjb),rhat);
    c1=corr([Fhatb(:,1) p1_cyt]);
    c2=corr([Fhatb(:,2) p2_cyt]);
    p1_cyb=sign(c1(1,2))*Fhatb(:,1);
    p2_cyb=sign(c2(1,2))*Fhatb(:,2);
    
    xb=[ones(rows(p1_cyb),1) p1_cyb p2_cyb cyb_ave1 cyb_ave2 cyb_ave3 cyb_ave4 cyb_ave5 cyb_ave6 dsb_1 inflb_1 inflb_2];
    [bb,bintb,rb,rintb,statsb] = regress(inflb_h,xb);
    Gb=nw((rb*ones(1,cols(xb))).*xb,nwl);
    Vb=inv(xb'*xb)*Gb*inv(xb'*xb);
    seb = sqrt(diag(Vb));
    betab(j,:)=(bb-b)';
    tstatb(j,:)=((bb-b)./seb)';
    
end;

fprintf(' h = %6.4f\n',h);
fprintf('\n estimation results for baseline inflation model\n')
fprintf('    coeff.     s.e.    boot s.e.  90%% bootstrap CI\n')
% equal-tailed percentile-t bootstrap
disp([b se std(betab)' b-se.*quantile(tstatb,0.95)' b-se.*quantile(tstatb,0.05)'])
% Hall's percentile bootstrap
%disp([b se std(betab)' b-quantile(betab,0.95)' b-quantile(betab,0.05)'])
fprintf(' Rbar^2 = %6.4f\n',R2);
%fprintf('\n correlation matrix of PCs and AVEs\n')
%disp(corr([p1_cy p2_cy cy_ave1 cy_ave2 cy_ave3 cy_ave4 cy_ave5 cy_ave6]))
fprintf(' R^2 = %6.4f\n',R2_2);

end